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Abstract 

We propose a formulation of adaptive computation of free energy differences, in the ABF or 
nonequilibrium metadynamics spirit, using conditional distributions of samples of configurations 
which evolve in time. This allows to present a truly unifying framework for these methods, and to 
prove convergence results for certain classes of algorithms. From a numerical viewpoint, a parallel 
implementation of these methods is very natural, the replicas interacting through the reconstructed 
free energy. We show how to improve this parallel implementation by resorting to some selection 
mechanism on the replicas. This is illustrated by computations on a model system of conformational 
changes. 
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I. INTRODUCTION 



One of the most important goals of molecular simulation is the computation of free energy 
differences as a function of some selected degrees of freedom, called reaction coordinates. 
The dynamics of the system can indeed often be split into slowly evolving degrees of freedom, 
which determine the reaction coordinates to be used, and other rapidly evolving degrees of 
freedom. The free energy differences allow to characterize global changes in the system 
under study, and give information about the relative stabilities of several species, as well as 
their transition kinetics. However, the free energy barriers to overcome are so large in many 
applications that a computation based on a straightforward sampling is unfeasible since the 
system remains stuck in metastable free energy sets. 

A classical technique to compute free energy differences is thermodynamic integration, 
dating back to Kirkwood^, which mimics the quasi-static evolution of a system as a succes- 
sion of equilibrium samplings, which amounts to an infinitely slow switching between the 
initial and final states. Another classical technique is the free energy perturbation method, 
introduced by Zwanzig^, which recasts free energy differences as a phase-space integral, so 
that usual sampling techniques can be employed. Notice also that there exist many refine- 
ments for those two classes of techniques, such as umbrella sampling^. 

More recently, methods relying on nonequilibrium dynamics have emerged. They follow 
the pioneering work of Jarzynski^, or use some adaptive dynamics such as the Wang-Landau 
approach^, the adaptive biasing force (ABF)^^, or the nonequilibrium metadynamics^. 
These approaches use the whole history of the exploration process to bias the current dy- 
namics in order to force the escape from metastable sets. This is done by simultaneously 
estimating the free energy from an evolving ensemble of configurations of the dynamics, and 
using this estimate to bias the dynamics, so that the effective free energy surface explored is 
flattened. In the long time limit, the bias exactly gives the actual free energy profile. Adap- 
tive methods could therefore be seen as umbrella sampling with an evolving potential. This 
was already noticed in a previous study presenting an adaptive dynamics as a 'self-healing 
umbrella sampling—. 

To present the adaptive methods mentioned above in a general and unifying framework, 
it is convenient, as is done in^, to consider ensemble of realizations (see Eq. (jl])). The 
system is then described by the distribution of the configurations of this ensemble in the 



2 



limit of an infinite number of replicas simulated in parallel. The key point is to reformulate 
the computation of the bias of adaptive dynamics, using conditional distributions (that is, 
distribution of the configurations for a given value of the reaction coordinate) of the latter 
sample. This was already proposed in^ in the equilibrium case, and is somewhat implicit in^. 
This concept clarifies the presentation of adaptive methods, allows mathematical proofs of 
convergence^^ or at least, existence of a stationary state of the dynamics (still in the case of 
an infinite number of replicas), and suggests natural numerical strategies: the discretization 
may be done through a parallel implementation of several replicas of the system, which all 
contribute to construct the free energy profile. Such a parallel implementation was already 
proposed in^^ in the case of metadynamics. We show here how an additional selection process 
on the replicas can enhance the sampling of the reaction coordinates in comparison with a 
straightforward parallel implementation. 

The paper is organized as follows. In Section [Tll we describe the general formalism for 
adaptive dynamics, using conditional probabilities, and show how to update the biasing 
potential in order to compute the free energy profile in the longtime limit, using a fixed- 
point strategy. Some applications of this formalism are presented in Section UTTt and allow to 
recover the usual adaptive dynamics such as the nonequilibrium metadynamics, the Wang- 
Landau scheme or the ABF method. We then discuss possible parallel implementation 
strategies in Section HVl In particular, it is shown in Section [I V Bl how a selection process can 
enhance the straightforward parallel implementation. This is finally illustrated by numerical 
results for a toy model of conformational changes in Section |Vl 

II. A GENERAL FRAMEWORK FOR ADAPTIVE DYNAMICS 
A. Notations 

For a system described by a potential V{q), the Boltzmann measure in the canonical en- 
semble is exp {—l3V{q)) dq (where Z is a normalization constant, the so-called partition 
function). Consider a reaction coordinate ^, taking values in the one dimensional torus, or 
in the interval [0, 1]. In the latter case, reflecting boundary conditions for the dynamics on 
the two extremal values C,{q) = 0, C,{q) = 1 are used. The free energy (or potential of mean 
force (PMF)) to be computed is defined up to an additive constant by the normalization 
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of a Boltzmann average of the configurations restricted to a given value of the reaction 
coordinate: 

A{z) = -p-' In j exp(-/?V(g)) S^^,)^,. (1) 

The PMF A exactly gives the Boltzmann weights of the equilibrium distribution of the 
reaction coordinate. The so-called mean force A'{z) along the reaction coordinate can also be 
expressed as the Boltzmann average of a real- valued "force" field over the configurations 
restricted to a given value of the reaction coordinate ^(g) = z: 



A{z) = . (2) 



Here and in the sequel, we denote by A' the derivative of A with respect to z. 

The force field can be expressed only in terms of derivatives of first and second order 
of the reaction coordinate ^ as 

For the mathematical derivation of this formula, and extension to multi-dimensional reaction 
coordinates or phase-spaces with holonomic constraints, we refer for example to^>^»^. 



B. Adaptive biasing dynamics 

Adaptive dynamics are defined through the dynamics used, which dictates the distribution 
of the configurations at equilibrium (see Section fllB ip . a biasing potential (Section 111 B2p . 
and the way this potential is updated (see Section 111 B 31 for a heuristic derivation in the 
equilibrium case motivating the general setting of Section IIIB4P . 



1. The dynamics 

Trajectories t ^ Qt are computed according to some dynamics which are ergodic with 
respect to the Boltzmann measure when the potential is time-independent. For instance, 
the Langevin dynamics or the overdamped Langevin dynamics may be used. We will de- 
note by ipti(l) the probability distribution (or density) of configurations at time t. This 
distribution will be used to update the biasing potential Abias- 
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From a practical point of view, when M replicas {Qt^)i=i,...,M of the system are simulated 
in parallel, the density of states iptiq) is approximated by the instantaneous distribution of 
the replicas 

1 

1=1 

In some cases, the density of states can also be approximated using the distribution of 
configurations along the trajectory, relying on some ergodic assumption. 

The definition of adaptive methods requires the definition of two important quantities 
obtained from the distribution iptio)- The first one is the distribution ip^ of the reaction 
coordinate values, which is the marginal law of ipt with respect to ^: 

V-t (^) = j S!^{q)-Z■ (5) 

This quantity will be useful to propose a biasing potential (see Eqs. ( !T2|) - (fT4|) ). Another 
important quantity is the conditional average of some function h for some fixed value of the 
reaction coordinate: 

Such averages are used to propose biasing forces (see Eqs. (fT^ - (IT^ ). 



2. The biasing potential 

In adaptive dynamics, the interaction potential is time-dependent: 

Vt{q) = V{q)-A^Ut,m)- (7) 

The biasing potential Abias? whose precise form varies according to the method under study, 
depends only on q through the reaction coordinate value ^{q) and is updated using the 
history of the configurations. It is expected that this biasing potential converges (up to an 
additive constant) toward the free energy A given by ([1]) in the long-time limit, so that the 
equilibrium distribution of the reaction coordinate is the uniform distribution. 

The key idea common to all adaptive methods is to resort to a fixed point strategy, in 
order for the observed free energy to converge to a constant or the mean force to vanish, 
and the dynamics to reach equilibrium (see the updates ([9]) or (fTT]) in the equilibrium case 
and f|T4|) or f|T5l) in the nonequilibrium case). 
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3. Updating the biasing potential - The equilibrium case 



To derive a possible form for the biasing potential, let us first assume that the system is 
instantaneously at equilibrium with respect to the biased potential Vt, i.e. Qt has density 
ipt'^iq) = Zf^ exp(— /5Vt(g)). In this case, resorting to (JT]), the observed free energy (see ( IT^ 
for a general definition) is 

- In j r^q) Sa<i)-z = A{z) - A^^t, z) + /^'^ In Z,. (8) 

Thus, for a characteristic time r to be chosen, an update of Abias of the form 

9* Aias(t, z) = -^\n f rt\q) Sa,)-z (9) 



is such that A'^[i^^{t) — > A' when t +oo exponentially fast with rate 1/r. Notice that 
we stated the convergence in terms of the mean force, because, in view of the constant 
term f3~^ iTiZt in Eq. ([8]), the potential of mean force only converges up to a constant to the 
true potential of mean force. 

Similar considerations hold for the mean force: replacing the potential V with Vj given 
by ©, and resorting to ©-([S]), the observed mean force (see f|T3l) for a general definition) 
is 

— = A'{z) — v4{,i3^g(t, z), (10) 
since /^'(g) = /^(g) - A'^,,,{t,^{q)). An update of A'^,,,it) of the form 

* ^"r /^r(g)%.)-.(rfg) ^ ^ 

is therefore such that A'^^i^^it) A' when t +oo exponentially fast with rate 1/r. 

4. Updating the biasing potential - The nonequililibrium case 

Now, in general, the system is not at equilibrium for the potential Vt'. ipt 7^ We 
use the above procedure as a guideline to update the biasing potential y4bias(^, z). To derive 
equations for the biasing potential, let us first define two quantities. The first one is the 
observed free energy or the observed potential of mean force, defined as 

^pot,obs(i, 2;) = -^"Mn j ^t{q)k{q)-z- (12) 



This quantity can be interpreted as the free energy associated with the ensemble of config- 
urations with density of states iptig) (see Eq. ([1])). The observed free energy v4pot,obs(^, ^) 
is high when the number of visited states with reaction coordinate value z is small. In the 
long-time limit, the distribution of the reaction coordinate is expected to be uniform, so 
that the observed free energy is constant. 

In the same way, the observed mean force is defined as the conditional average of the 
time-dependent biasing force for a given value of the reaction coordinate: 

A> (f .^- I ^^'(^) ^*(^) ^^fa)-- - / ^^(g) '^^W-^ A' (, ,\ 

This quantity can be interpreted as the mean force associated with iptil) (see Eqs. ([2])- ([3])), 
minus the biasing force at time t. It is expected to vanish in the long-time limit, so that the 
corresponding observed free energy is also constant. 

The fixed point strategy relies on two different ways of updating the bias (the updating 
functions Ft and Gt are increasing functions such that Gt{0) = 0): 

• The first strategy, which may be called Adaptive Biasing Potential (ABP) method, is 
inspired by ([9]). The bias is updated in its potential form, preferably increased (resp. 
decreased) for reaction coordinate values such that the observed free energy is high 
(resp. low): 

(ABP) 9jAbias(t,2:) =Fj(Apot,obs(t,;2)); (14) 

• The second strategy, the usual ABF method, is inspired by ffTTj) . The bias is updated 
through the mean force: the biasing force is increased (resp. decreased) for reaction 
coordinate values such that the observed mean force is positive (resp. negative): 

(ABF) dtA'^;^{t, z) = GMLcco^sit, z)). (15) 

Let us emphasize at this point that the ABF and the ABP methods yield very different 
biasing dynamics, since the derivative of ( fT2l) with respect to z is different from ( fT3l) (This 
would not be the case if the system was at equilibrium: the derivative of ([9]) with respect to z 
is equal to ( ITT]) ). This difference becomes critical for multi-dimensional reaction coordinates, 
where the biasing force no longer derives from a potential in general. 
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C. Consistency of the method 

Let us show that within this formahsm, any stationary state of the ABP or ABF methods 
gives the true mean force A' to be computed (and therefore the true PMF up to an additive 
constant). Recall at this stage that we are dealing with distribution of replicas, which 
arise formally in the limit of an infinite number of replicas. For a stationary state where 
the biasing potential has converged to y4bias(oo), the ergodicity property of the dynamics 
ensures that samples of configurations of the system are distributed according to ipoo = 

Z-^exp[-P{V-AuUoo,m- 

The observed free energy or mean force given by Eqs. 0121) and fll3p then both verify 

^pot,obs(oO'2;) = A'f^^^^^^^^{oo,z) = A'{z) - A'^^^{oo,z). The updating equations Eqs. ^ 

and (|T5|) yield respectively 

F^{A{z)-AuUoo,z)) = 0, (16) 

Goo(A'(2)-A;,,,(oo,^)) = 0, (17) 

so that (taking the derivative with respect to z in (fT6|) ). A[^^^{oo) = A' in both cases thanks 
to the strict monotonicity of the updating functions. Let us also notice that, at convergence, 
the values of the reaction coordinate are distributed uniformly: J ipooil) ^^{q)-z = 1- 

However, let us emphasize that we did not give any convergence result at this point. We 
merely showed that, if the dynamics converges, then the limiting state is the correct one. 
To prove convergence starting from an arbitrary initial distribution is a difficult task, and 
can only be done for certain dynamics (see the corresponding results in Section [ITI|) . 

III. APPLICATION TO USUAL ADAPTIVE DYNAMICS AND CONVER- 
GENCE RESULTS 

We present in this section some applications of the formalism of Section [Tll and show 
that the usual adaptive methods can indeed be recovered. This is summarized in Table [U 
which gives a classification of adaptive methods. We then give a rigorous convergence result 
for some class of adaptive methods. 

Table I 
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A. Met adynamics 



Adaptive strategies can be used with met adynamics^. The configuration space is ex- 
tended by considering an additional variable z representing the reaction coordinate, and 
the dynamics is denoted t ^ {Qt,Zt). The associated extended potential incorporates a 
coupling between this new variable and the reaction coordinate ^: 

V>^{q,z) = V{q) + f^{z-aq)f, 

for some (large) /i > 0. In this case, the new reaction coordinate considered is ^meta('?5 z) = z 
and the free energy is thus given by: 



A^{z) 



-/?"Mn j exp{-pV>'{q,z))dq. 



It is easy to check that, up to an additive constant, A'^ ^ A as /i — > +00, with A given 
by ([I]). The adaptive strategies presented above applied to this extended dynamics allow 
to recover the free energy A^. The corresponding dynamics may be called meta- Adaptive 
Biasing Potential (m-ABP) and meta-Adaptive Biasing Force (m-ABF) methods. 

Strategies relying on biasing potentials are reminiscent of fiooding strategies^ such as 
the nonequilibrium metadynamics^. An example of an m-ABP method is metadynamics^^ 
when the biasing potential is applied to the extended variable. The updating function 
does not depend on time and is given by Ft{x) = — 7exp(— /5x) for some constant 7 > 0. 
The ensemble of configuration used in the adaptive update is obtained from M replicas 
{Qt^\ Zl'^) running in parallel, so that 

1 

The resulting biasing potential at time t penalizes the values of the reaction coordinate 
already visited according to (see (fl^ ): 

AHas(t, Z) ^ <L(t, ^) = -Zj2 ^Z^'-z (18) 

j=l 

In the case of an overdamped Langevin dynamics with M = 1 for example, the resulting 
equations of motion are therefore: 

dQt = -VV{Qt) dt + ^l{Zt - e(Qi))Ve(gt) dt + ./W^dW?, 
dZt = -fi{Zt - e(Qt)) dt + y^dWf - 7V, ( £ 6z,-z ds^ dt, 



where the processes W^, Wf' are independent standard Brownian motions. The usual 
met adynamics are recovered when, as in^ii^, the Dirac masses Szt-z are discretized in the 
last equation and in ffTSl) using Gaussian functions. We also refer to^ for an error analysis. 



B. The Wang-Landau algorithm 

Another famous instance of an ABP dynamics, usually defined in discrete spaces, is the 
Wang-Landau algorithnt^. The biasing potential is constructed in a similar fashion to f|T8l) . 
without extending the configuration space and with only one replica. The updating function 
is modified during time as F^i^x) = — 7(t) exp(— /5a;), so that 

Ahiasii, = - l{s)SaQs)-zds. (19) 
Jo 

If 7(t) -^0 slowly enough, it is possible to prove the convergence of the dynamics, the rate 
of convergence of 7(t) being controlled by the nonuniformity of the histogram of the time 
distribution of the reaction coordinate (see^ for more precisions on the convergence results) . 



C. The ABF method 

The usual ABF bias^ is given by averaging the local force over the configurations 
visited by the system. It is recovered in the formalism we propose by considering one replica 
of the system, and an updating function of the form Gt{x) = 7X in the limit 7 oo. This 
gives indeed: 

Since there is only one replica, the density ipt{s) is approximated by a trajectorial distribu- 
tion, for example 

MQ)-7f / ^Qs-.ds (21) 
for some averaging time T > and t > T. 



D. A rigorous convergence result in the ABF case 

A rigorous convergence result of the ABF algorithm with the update (l20l) can be shown 
in the case of an overdamped Langevin dynamics^. It applies in the limit of an infinite 
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number of replicas simulated in parallel using (jl]). Let us detail this point. Consider the 
modified overdamped Langevin dynamics: 

dQt = -V{V + 2p-Hn |Ve| - Abias(t,0)(Qt) \V^\-\Qt)dt+^/2^\V^r\Qt)dWt,{22) 

with the update ([20D: A'^i^^{t,z) = {f^)t,z- The process Wt is the standard Brownian 
motion. This dynamics is the usual overdamped Langevin dynamics for the potential Vt when 
|V^| = 1. Notice that in the case of a metadynamics-like implementation ('m-ABF'), the 
modified dynamics is actually the usual overdamped Langevin dynamics since ^meta(Q') ^) — ^ 

and thus |V^meta| = 1- 

The reason we propose to add the three terms depending on |V.^| is twofold: for the 
dynamics (l22!) - (!20l) . it can be proved^ that (i) the distribution ip^ of the reaction coordinate 
satisfies 

which is a pure diffusive behavior for the marginal law of the reaction coordinate. In par- 
ticular, the initial metastable features of the free energy landscape are overcome in a time 
scaling as the square of the length of the reaction coordinate set; (ii) the observed mean 
force ^bias(^) converges to A' when t +00. The convergence is exponential, the rate of 
convergence being limited by the minimum between the rate of convergence in each subman- 
ifold ^{q) = z (see^ for more details), and the rate of convergence in the reaction coordinate 
space. Let us emphasize that this proof of convergence does not assume that the system is 
at equilibrium at any time. 

IV. PRACTICAL IMPLEMENTATION STRATEGIES 

Relying on the definition (jlj) of the distribution of configurations, adaptive dynamics 
can be easily parallized by using a large number M of replicas that interact through the 
biasing potential or the biasing force. We first show in this section how to discretize the 
dynamics and the biasing potential (section llV Ap . and then, how this implementation can 
be improved using some selection process (section HVB]). 
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A. Discretization of the biasing potential 



In order to compute in practice the conditional or marginal distributions needed to update 
the biasing potential, there are basically two approaches, relying either on ergodic limits or 
on ensemble averages. Both approaches may be combined in practice in order to obtain 
smooth profiles. For example, when only a limited number of replicas M is used, the density 
ipt{<l) given by (jlj) is not regular, and some local averaging is necessary (see e.g. Eq. fl23l) or 
Eq. dH). 

We detail the implementation in the ABF case for example. The ABP case can be treated 
in a similar way (see also^^). The instantaneous conditional average of some function h is 
typically approximated by 



M 
t,z 



where Qt'^ is the i-th replica at time t and 61 is some approximation of the Dirac distribution 
6z, such as a gaussian function with standard deviation e or the indicator function of an 
interval of size e. In order to regularize these averages over the replicas, some time averagings 
may be used (as in ([21])) such as 



ds 



t.z 



ds 



or 



Kr{t 



ds, 



(23) 



(24) 



with a convolution kernel KT-(t). For instance, KT-{t) = 1j>ot ^e Many other regular 
izations relying on a (local) ergodicity property could of course be used. 



B. Enhancing the sampling through a selection process 

A general strategy to improve the straightforward parallel implementation (jl]) is to add 
a selection step to duplicate "innovating" replicas (replicas located in regions where the 
sampling of the reaction coordinate is not sufficient), and kill "redundant" ones. One way 
to perform an efficient selection is to consider an additional jump process quantified by 
a field S{t, z) over the reaction coordinate values. Each replica trajectory (Q!.'*^) is then 
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weighted by exp{J^ S{s,^{Ql'^)) ds), which naturally gives birth/death probabilities for the 
selection mechanism, in the spirit of Sequential Monte Carlo (SMC) methods^S, or Quantum 
Monte Carlo methods (QMC)^. A possible choice is 

S = c^, (25) 

where c is a positive constant. This method thus enhances replicas in the convex areas of 
the density tp^, where free energy barriers still need to be overcome. When convergence is 
reached, ip^ is uniform and the selection mechanism vanishes. 

When the selection step is used with the overdamped Langevin dynamics (l22l) . it can 
be shown that the distribution of the reaction coordinate values still satisfies a simple 
diffusion equation, but with a higher diffusion constant: 

This method thus enhances the diffusion in the reaction coordinate space, but the conver- 
gence rate is still limited by the relaxation in each submanifold ^(g) = z. 

The jump process can be computed in practice by attaching a birth time and a death 
time to each replica. The death time is decreased when the source term S given by ( l25l) is 
positive, otherwise the birth time is decreased. When the death time is zero, the replica is 
replaced by another replica chosen at random. The birth process is handled in a similar way: 
when the birth time is zero, a replica chosen at random is replaced by the replica whose birth 
time is zero. We refer for example tc3 for more precisions on the practical implementation 
of such a procedure, as well as other strategies to handle birth/death processes. 



V. NUMERICAL RESULTS 



We finally present an application of the selection strategy proposed in section IIVBI to a 
model system of conformational change in solution. We consider a two-dimensional system 
composed of particles in a periodic box of side length I, interacting through the purely 
repulsive WCA pair potential^'^: 



VwcA(r) 



4e 
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(7\6 

r 



+ e if r < To, 
if r > tq. 
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In these expressions, r denotes the distance between two particles, e and a are two positive 
parameters and Among these particles, two are designated to form a solute dimer 

while the others are solvent particles. For these two particles, the above WCA potential is 
replaced by a double-well potential 

(r — ro — w)' 



V.(r) = h 1 



where h and w are two positive parameters. The potential Vs exhibits two energy minima, 
one corresponding to the compact state r = tq, and one corresponding to the stretched state 
r = ro + 2w. The energy barrier separating both states is h. The reaction coordinate ^ is 
the bond length r of the dimer. 

In practice, the Dirac distribution are approximated by indicator functions of intervals 
of size Az = 0.05. The parameters used for these computations are = 16 particles, at 
particle density p = N/t^ = 0.25a-'^, a = 1, w = 0.7, e = 1, /i = 20 and /3 = 5. We 
consider M = 2000 replicas evolving according to an overdamped Langevin dynamics, with 
a time step At = 10~^. The reference computation is done with M = 5000 replicas and the 
reference mean force profile is obtained by averaging the profiles on the time interval [5, 10]. 
The profiles are regularized in time by using ( l24l) with r/At = 100. The initial conditions 
are such that the dimer bond lengths of all replicas are close to tq. We consider in the 
sequel the interval [20,^1] = [1.1,2.55] (since, with the parameters chosen here, ro — 1.122, 
ro + 2w ~ 2.522 and Az = 0.05), containing n = 30 bins. 

We present in Figure [1] free energy difference profiles (averaged over K = 100 indepen- 
dent realizations) obtained with the parallel ABF dynamics ( l20l) . with and without the 
birth/death selection term (l25l) (with c = 10), at a fixed time tfiguro = 0.1. The standard 
deviation of the profiles {A[, . . . , A'j^) for K independent realizations is 



\ 



k=l 



where A'{z) = Ylk=i mean force averaged over all the realizations. The 

associated 95% confidence intervals (or errors bars) are 



1 96 1 96 

A'{z) - -j=aA'{z). Al{z) + -j=OA'{z) 
y K y K 



(26) 



The curves plotted in solid lines in Figured] are the averages A!, and the curves plotted in 
dashed lines are A!_ and A!j^. Notice that the mean force profile obtained when the selection 

14 



process is turned on is converged (since the curves A', A'_, A'^ and the reference curve are 
almost indistinguishable) . 

Figure I 

The comparison shows that the selection process improves the rate of convergence of the 
algorithm and accelerates the exploration process on the free energy surface. Indeed, the 
profile obtained when the selection process is turned on is quickly really close to the reference 
profile. On the other hand, with a straightforward parallelization, only a small fraction of 
replicas has escaped from the initial free energy metastable state at time tfigurc to explore 
the free energy metastable set corresponding to bond lengths around tq + 2w. 

To precise these qualitative features, we further perform two quantitative studies for 
several values of c: 

(i) Table [TTl and Table IIIII make precise the convergence of the profiles to the reference 
profile in a quantitative way. The measure of error we consider is 

SA= max \A{z) — Aj.c{{z)\, 

Z()<Z<Zl 

where Aj-ef is the reference profile, and A{z) = f^^ A' is the averaged potential of mean 
force, obtained as the integral of the mean force averaged over all the realizations. In 
practice, we consider the following approximated deviation between PMF profiles: 



6 An = max 

0<i<n 



Az. (27) 
A 95% confidence interval is obtained as with 

Az. 



5^ An = max 

0<j<n. 



j=i V 



(ii) Figure [2] presents the fraction of replicas which have crossed the free-energy barrier 
(averaged over the K = 100 realizations), i.e. the instantaneous fraction of particles 
such that r > vq + w. Notice that we expect this fraction to converge to 0.5 (up to 
some errors due to statistical fluctuations and to the binning of [2:0, ^i])- 

Table II 
Table III 
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Figure II 

As can be seen from the different escaping profiles of Figure [2], the selection process really 
accelerates the transition from one free energy metastable state to the other. This is due 
to the fact that the birth and death jump process triggers non local moves, as opposed to 
the traditional diffusive exploration of adaptive dynamics. The numerical results of Table [ITl 
show that it is very interesting to consider a selection process, especially at the early stages 
of the simulation. This selection is even more efficient when the number of replicas increases 
(see Table [In|) . In conclusion, the selection process seems to be an efficient tool to improve 
the exploration power of the adaptive dynamics. 
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Table and Figures captions 



• Table I. Classification of adaptive methods. 

• Figure I. (color online) Free energy difference profiles obtained with the parallel ABF 
algorithm (in reduced units), for a time tfigure = 0.1 and averaged over K = 100 
independent realizations: with birth/death process (c = 10, red), without birth/death 
process (blue), reference computation (black). Solid line: average mean force; dashed 
lines: upper and lower bounds of the 95% confidence intervals (see Eq. (12^1)). 

• Table II. Deviation SAn from the reference PMF profile (given by Eq. ( !27ll ) as a 
function of the selection parameter c (c = when the selection is turned off) and the 
simulation time tgimu- The 95% confidence interval is given in brackets 
{K = 100 realizations). 

• Table III. Deviation 6An from the reference PMF profile (and associated error bars) 
when c = 10 for different number of replicas {K = 50 realizations). 

• Figure II. (color online) Average fraction of the replicas in the region r > vq + w a.s 
a function of time, for c = (no selection, black), c = 2 (blue), c = 5 (green), c = 10 
(red). 
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Adaptive Biasing 


Adaptive Biasing 






Force (OtA'^.J 


Potential (^t^bias) 




Dimension n (V) 




ABP^ 


Dimension n + 1 (V^'^) 


m-ABF 


m-ABP^>i^ 


TABLE I: Lelievre et al., Journal of Chemical Physics. 


c 


^simu — 0.05 


0.1 0.2 


0.4 





9.51 (7.73-11.3) 


18.0 (14.8-21.2) 19.5 (18.3-20.7) 


0.066 (0.056-0.075) 


2 


20.4 (17.0-23.8) 


5.69 (5.55-5.82) 0.020 (0.016-0.023) 


0.034 (0.029-0.038) 


5 


22.9 (20.9-24.9) 


0.22 (0.19-0.25) 0.027 (0.022(0.032) 


0.026 (0.022-0.031) 


10 


10.4 (10.4-10.4) 


0.035 (0.029-0.041) 0.028 (0.023-0.032) 


0.032 (0.027-0.037) 


TABLE II: Lelievre et al., Journal of Chemical Physics. 


number of replicas 


tsimu = 0.05 0.1 


0.4 




1000 


23.3 (20.4-26.3) 0.45 (0.39-0.50) 


0.064 (0.054-0.074) 




2000 


11.2 (11.2-11.2) 0.034 (0.025-0.042) 


0.032 (0.024-0.039) 




10000 


2.05 (1.54-2.56) 0.026 (0.019-0.033) 


0.022 (0.016-0.028) 



TABLE III: Lelievre et al., Journal of Chemical Physics. 
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FIG. 1: Lelievre et al., Journal of Chemical Physics. 
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FIG. 2: Lelievre et al., Journal of Chemical Physics. 



